Confinement modifies interactions of ultracold dipolar gases on optical lattices 
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We study the effective dipole-dipole interactions in ultracold quantum gases on optical lattices 
as a function of imbalance in confinement along the principal axes of the lattice. We demonstrate 
that the effective dipole-dipole interaction in the lattice decays exponentially with the inter-particle 
separation at short to medium distance on the lattice scale and has a long-range power-law tail, 
in contrast to the pure power-law behavior of the dipole-dipole interaction in free space. The 
effect can be sizable; we identify differences of up to 36% from the free-space interaction at the 
nearest-neighbor distance in quasi-lD arrangements. The modified interaction arises from quantum 
fluctuations induced by heavy tails of the localized single-particle probability distributions, and 
also relies crucially on imbalance in confinement, due to the d-wave anisotropy of the dipole-dipole 
interaction. Using matrix product state simulations, we demonstrate that use of the correct lattice 
dipolar interaction leads to significant deviations from many-body predictions using the free-space 
interaction in the lattice. Our results are relevant to up and coming experiments with heteronuclear 
molecules, Rydberg atoms, and strongly magnetic atoms in optical lattices. 



Recent experimental progress in cooling heteronuclear 
polar molecules with large electric dipole moments [T], 
Rydberg atoms [2] , and atoms with large magnetic dipole 
moments, in particular Chromium [3], Erbium [4], and 
Dysprosium [5] , has sparked interest in the properties of 
ultracold dipolar gases. While in many ultracold atomic 
systems interactions are short-range and well-modeled 
by a contact pseudopotential, the interactions in dipolar 
gases have a long-range and anisotropic character, de- 
caying as 1/r 3 with the separation r between particles. 
These features of the dipole-dipole interaction have lead 
to a variety of intriguing theoretical proposals such as ex- 
otic pairing and bound states in ladder geometries [SJ [7] 
and the realization of quantum liquid crystal states of 
matter [9] . Even for atoms with relatively weak dipole 
moments, such as Rubidium, dipole-dipole interactions 
can play a significant role [10] . In this Letter, we show 
that the effective dipolar interaction in a lattice is not 
actually 1/r 3 as commonly believed, but is in fact ex- 
ponential at moderate separations, and only behaves as 
1/r 3 for large separations. 

A key component of our analysis is the presence of a 
continuous, periodic potential, which for ultracold atomic 
and molecular gases is provided by an optical lattice. As 
first discovered by Kohn [11], the most localized set of 
orthogonal single-particle states with the symmetries of 
the lattice generally feature an exponential decay. We 
find that the effective dipolar interaction in the lattice 
depends essentially on the exponential tails of the single- 
particle states rather than only on their widths. Hence, 
approximating the localized states with functions which 
match only the mean width will fail to accurately capture 
the effective interaction in the lattice. We will show that 
corrections of order 36-48% arise for interactions at the 
nearest-neighbor distance in moderately confined quasi- 
low-dimcnsional scenarios. 

An additional essential ingredient for our findings is an 
imbalance in the degree of confinement along the princi- 



pal axes of the lattice due to the anisotropic character 
of the dipole-dipole interaction. In this work, we char- 
acterize confinement using the curvature of a lattice site 
minimum. Our work builds on a wealth of confinement- 
induced phenomena in ultracold quantum gases, such as 
confinement-induced resonances |12) . the fermionization 
of a ID Bose gas [T3], and the Berezinskii-Kosterlitz- 
Thouless transition in a quasi-2D Bose gas [2] . In dipo- 
lar gases, the effects of confinement have been studied 
in harmonic traps |15[ I16j . There, it has been shown 
that strong confinement along the axis of a field orient- 
ing the dipoles leads to purely repulsive interactions in 
the weakly confined plane. Additionally, because of the 
anisotropic character of the dipole-dipole interaction, the 
stability of a dipolar BEC displays a strong dependence 
on anisotropy in external confinement, and stable solu- 
tions can even take surprising forms such as the "red 
blood cell" dipolar BEC [17] • Similarly, anisotropic lat- 
tice confinement has surprising differences from precon- 
ceptions based on uniform isotropic systems. Dipolar 
interactions in confined geometries appear in many other 
branches of physics, such as in ferromagnetic nanostruc- 
tures [IB] , where the dipole-dipole interaction plays a key 
role in the dispersion relation for spin waves in ferromag- 
netic films pjj] and wires [20] . 

An important application of our results is in deriving 
many-body models to describe the low-energy physics of 
dipolar gases on lattices. Previous derivations, for exam- 
ple Refs. [2TM24) . assume that the interaction between 
localized states in the lattice has the same functional 
form as in the continuum. Since this amounts to replac- 
ing the localized single-particle probability distributions 
with delta functions, we will call this the delta function 
approximation (DFA). Performing matrix product state 
simulations [25 on infinite one-dimensional (ID) lattices, 
we demonstrate that the DFA can lead to significant er- 
rors in the determination of the phase diagram. 

For particles subject to a periodic lattice potential, the 
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energy eigenfunctions are Bloch functions ^nk(r) charac- 
terized by a quasimomentum index k in the first Brillouin 
zone (BZ) and a band index n. A more appropriate ba- 
sis for describing strong correlations in lattices is pro- 
vided by Wannier functions, which are the quasimomen- 
tum Fourier transforms of the Bloch functions, 

wt(r) = w(r - Ti ) = L- 1 / 2 £ keBZ e^'Mr) • (1) 

Here L is the total number of unit cells in a lattice with 
periodic boundary conditions and denotes the coordi- 
nate of lattice site i. For simplicity, we have restricted 
our attention to the lowest band. The extension of our 
results to multi-band situations is straightforward. 

It was shown in a seminal work by Kohn that the 
phases on the Bloch functions can be chosen such that the 
associated Wannier functions are exponentially decaying 
away from zero for one-dimensional centro-symmetric lat- 
tice potentials. This result has been extended to general 
one-dimensional lattice potentials [26], non- degenerate 
bands in arbitrary dimensions [2 7) , and general two- and 
three-dimensional insulators with vanishing Chern num- 
ber [55] . Hence, the exponential localization of Wannier 
functions, which will play an important role in our re- 
sults, is a general property which does not require fine 
tuning or a specific lattice structure. The numerical re- 
sults presented in this work employ a simple cubic lattice 
potential commonly used for ultracold quantum gases, 

V(v) = Y, ve{x , y , z} V u sm 2 (-Kv/a), (2) 

where a is the lattice constant. 

With the identification of the Wannier functions as the 
appropriate single-particle basis for describing strongly 
interacting particles in a lattice, we derive a many-body 
lattice model using the well-known procedure of expand- 
ing the field operator ip(r) in the complete basis of Wan- 
nier functions and substituting this expansion into the 
second-quantized expression for the interaction Hamilto- 
nian, 

H mt = ydrJdr'ft(v)i,^r')V int (r-r')i,(v')$(v), 

(3) 

where V[ n t is the two-particle interaction potential. The 
expansion of Eq. ^ in lowest band Wannier functions 
yields 

Hint = I y],- i a l &t a-.ia-j , (4) 

iiil 2, ^^lil2l 2 li 1 1 1 2 1 2 1 1 11 12 L 2 l l ' v ' 

where <x\ destroys a particle in a Wannier state centered 
at site i and 

W, ll3lili = /dr/ < ^/ ll , 1 (r)V irt (r-r')/ Ia ,,(r). (5) 

Here, /ii'(r) = w\{v)w\i(v) is a product of Wannier func- 
tions. The matrix elements U^i^v , which we will call 
Hubbard parameters, describe the interactions between 
particles localized in Wannier states. 



It is convenient to write the dipole-dipole interaction 
potential as the contraction of two rank-two spherical 
tensors [29] . 

V dd (K) = -# E' = - 2 (-l) 9 ^ 2) (R)[di ® d 2 fl , (6) 

(2) 

where R is relative position of the two dipoles, C q (R) 
is an unnormalized spherical harmonic, 

[di ® d 2 ]< 2 > = £ m C^^ 2 _ m;9 (d 1 ) m (d 2 ) ? _ m , (7) 

^rnimL-m = (.h m ih m 2\j'm) is a Clebsch-Gordan coeffi- 
cient, and (di) m represents the projection of the dipole 
operator for the i th particle along space-fixed spherical 
direction e m with e the z axis. For simplicity we as- 
sume that the expected dipole moments are position- 
independent. In this case, the expansion Eq. Q reads 

v v 2 j-\yg AA :\., [d x ® d 2 ] (2 ]a! a\ a v a v , 

4 ^lil2l 2 li *—iq=— 2 \ I ^liial^l^ l 1 ZJ —q 11 12 1 2 1 1 ' 

where the geometric contribution to the Hubbard param- 
eter is given by the integral 

offiV =- 2 /*/*'W')fe 2 W'')' ( g ) 

We use two different techniques to evaluate Eq. ([8| nu- 
merically. The methods and their convergence behavior 
are described in the supplemental material [3T]. We es- 
timate conservatively at most a 1% relative error in our 
presented results. 

For ultracold molecules in a DC electric field oriented 
along eo, only the q — terms in Eq. ^ are relevant 
in typical cases due to large energetic barriers to other 
dipole- allowed processes. Hence, while our methods ap- 
ply for all components of the dipole-dipole interaction, 
we focus on the q — components. Taking into account 
the statistics of the particles and neglecting assisted tun- 
neling terms [32], the expansion Eq. Q may be written 
for the dipole-dipole interaction as 

flint = U[~ ^/ «i("i - 1) + 2 ^h'shihi,] , 
i i#i' 

where a 3 U = [di <g> d 2 ][, 2) y/3/2, 

j 3^dd;0 t 3iv>dd:Q . ^ddrOi / n \ 

h = a Goobo > J i'A = a [^ii'i'i ± ^ii'ii'] ■ ( 9 ) 

The plus (minus) sign on the exchange term in Eq. ^ 
refers to bosons (fermions). 

With respect to the length scale of the lattice constant, 
the short range physics is given by the interactions la 
which occur within a particular unit cell. A point of com- 
parison for the on-site dipolar interaction Iq is provided 
by the dimensionless interaction energy Jq° = (H dd )/U 
of two particles in the ground state of a harmonic trap. 
For consistency, we choose that the oscillator frequency 
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matches the local curvature of a lattice minimum, giv- 
ing the condition l v = a/(irVv ). Here, l v is the har- 
monic oscillator length along Cartesian direction v and 
V v = Vu/Er with E R — h 2 n 2 /2ma 2 the recoil energy. In 
the case where the oscillator lengths in the xy plane are 
equal, t x = £ y = i±_ , we have that 



i 3 jho 



\/2[2/3 + (i 2 - /?(! - a 2 )- 1 cot" 1 ^)]/^ , (10) 



where £ = (izi 2 ^) 1 ^ 3 /a is the geometric mean oscillator 
length in units of a, a = £ z /£± measures the confinement 
anisotropy, and /3 = a(l — a 2 ) -1 / 2 . Notably, the interac- 
tion energy Jj° vanishes for isotropic confinement, a = 1. 
For a < 1, corresponding to stronger confinement along 
the quantization axis, the negative part of the dipole- 
dipole potential is suppressed. Hence, Jj° is positive for 
a < 1. In contrast, for a > 1 where confinement is weak- 
est along the quantization axis, Iq° is negative. 

For the effective interaction in a lattice potential, we 
first consider the case where all three lattice heights in 
Eq. Q are equal, V = V x = % = V z . We find that 
the on-site dipolar interaction Iq vanishes for any value 
of V > 2, which parallels the result for the harmonic 
oscillator. For the moderate- to long-range interactions 
Ii< t i, let us define Ij = A+jj with j along the x direction. 
The nearest-neighbor interactions I\ are within 1% of 
the DFA prediction j~ 3 for lattices of depth V > 7. For 
shallower lattices, the effective interaction is less, with 
the nearest-neighbor interaction I\ w 0.97 at V = 4 and 
ii s» 0.79 at V = 2 for bosons. In this weak-lattice 
regime, the exchange term in Eq. ^ plays a significant 
role, while it is essentially negligible for V > 7. 

To demonstrate the effects of imbalanced lattice con- 
finement, we evaluate the parameters I- } for two quasi- 
low dimensional scenarios. In the quasi-2D scenario, we 
take the z direction to be tightly confined with a lattice 
strength V z = 45 and the lattice strengths along the x 
and y directions to be equal, V± = V x = V y . In the 
quasi- ID scenario we take both the z and y lattices to 
be tightly confining with V z = V y = 45, and call the x 
lattice strength Vj_. These quasi-dimensional reductions 
are quite moderate and are commonplace in current ex- 
periments [3D]. In the confined scenarios, we find correc- 
tions to the DFA with larger deviations for smaller j and 
smaller V±. Here and throughout the rest of this paper, 
we measure j along the x direction. We characterize the 
modified dipolar interaction by fitting the numerically 
obtained data for Ij to the form 



Ij = a e exp(-b e j) 



wj 



(11) 



for j € [1,7], i.e., out to seven sites. The DFA pre- 
dicts no exponential, a e = or b e — > oo, and w = 1 
and p = 3 for the long-range contribution. We chose the 



ansatz Eq. (Ill from several fit functions because it had 
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FIG. 1: (Color online) Confinement-induced modifications of 
dipole-dipole interactions in quasi-lD and quasi-2D. In all 
panels the solid (dashed) lines denote the best fit parameters 
to Eq. (Ill in the quasi-lD (quasi-2D) scenarios as a function 
of the quasi-low dimensional lattice height V±. Fit parame- 
ters for (a)-(b) bosons and (c)-(d) fermions. Top panels are 
exponential weight a e (red) and decay constant b e (green). 
Bottom panels are percent differences of long-range weight w 
(red) and power p (green) with respect to the DFA. 

length scale a c ~ a/b e of the confinement-induced mod- 
ifications. However, we do not propose that the form 
Eq. (11) is exact. Across a wide range of quasi- low di- 
mensional confinement, a c ~ 0.2a, and so the moderate 
range over which confinement modifies interactions is a 
few lattice sites. 

The fit parameters for our numerically generated data 
are shown in Fig. [l] Here the solid (dashed) lines refer 
to the quasi-lD (quasi-2D) scenarios. The top panels are 
the short-range parameters a e and b e in Eq. ( |11[ ), with 
panel (a) (panel (c)) pertaining to bosons (fermions). 
The bottom panels are the percent differences of the long- 



range parameters w and p in Eq. ( 11 ) with respect to the 



DFA predictions. Here, panel (b) (panel (d)) pertains to 
bosons (fermions). As opposed to the V x = V y = V z case, 
we find that when V± < V z the effective interaction Ij is 
always greater than predicted by the DFA. Additionally, 
Iq is positive for V± < V z and is smaller than the inter- 



file lowest fitting error and it provides a characteristic 



action energy of the harmonic oscillator, Eq. (10) 31J. 
Finally, we note that for quasi-low dimensional confine- 
ment V± > 7, the predictions for bosons and fermions 
are the same to a few percent. This implies that the 
exchange contribution in Eq. Q plays no role for deep 
lattices. 

A semiclassical explanation for the increased effective 
interaction at small separation j is that a greater range 
of solid angle is included in the integration Eq. ^ due to 
the finite width of the function fu(r). However, if Eq. ^ 
is computed replacing the distribution fa(r) by a normal- 
ized rectangular distribution with anisotropic widths, the 
increase in Ij is drastically underestimated, often by two 
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orders of magnitude, even when the width of the distri- 
bution is much greater than the mean width of /u(r). 
Hence, we attribute the increased interaction not to the 
width of /ii(r), but to its exponential tails. Stated differ- 
ently, the increase in interactions is due to large quantum 
fluctuations from the heavy-tailed probability distribu- 
tion /ii(r). Additionally, as the dipole-dipole interaction 
selects out the part of the two-particle wave function with 
relative d-wave symmetry, the anisotropy in confinement 
is also vital to the modification of the interaction. 

To illustrate the implications of our findings for many- 
body physics, we study a model of quasi-lD hard-core 
bosons with long-range dipolar interactions |33| 



J7 = -i£ 



(ij 



■ 'j *3 



(12) 



Here, the nearest-neighbor tunneling amplitude is t, (i, j) 
denotes nearest- neighbor pairs i and j, Ij-i = Ij t i is the 
dipole-dipole interaction for a separation of (j — i) sites, 
and p is the chemical potential. We compute the phase 
diagram of Eq. (12 1 for two realizations of Ij using the 



infinite size variational ground state search algorithm for 
matrix product states (iMPS) |34j . In one realization, we 
use the DFA, Ij = j~ 3 . In the second realization, Ij is the 
actual effective interaction in the lattice, Eq. Long- 
range interactions are facilitated within iMPS by using 
matrix product operators and fitting the interaction Ij 
to a sum of exponentials [35] . We take the lattice heights 
to be V z = V y — 45, V x = 6, which fixes the tunneling en- 
ergy t. The coefficient U can be tuned for ultracold polar 
molecules by an applied DC electric field. One can deter- 
mine the strength of the interaction U knowing only the 
expected dipole moment and the lattice constant from 
single-particle physics, and so scaling the phase diagram 
to U is appropriate. 

The phase diagrams of Eq. (12 1 for the DFA and 



lattice-modified dipolar interactions are shown in Fig. [2] 
The bordered areas represent predicted regions of crys- 
talline phase (CP) with density p = 1/2 and a non- 
vanishing single-particle gap. The remainder of the plot 
is a gapless superfluid (SF) phase. In the CP, the density 
correlation function J\f(r) = ((jiq — p)(n r — p)) behaves as 
Af(r) —> cst.(— l) r as r — > oo and the single-particle den- 
sity matrix A(r) = (aja r ) is exponentially decaying with 
r. The area bounded with dashed lines represents the 
region of CP computed using the DFA in the rescaled pa- 
rameters p. — p/U, t = p/U. The hatched area bordered 
with solid lines represents the CP boundaries computed 
with the actual lattice dipolar interaction and p = p/U, 
t = p/U. The region predicted by the DFA is shifted 
both in chemical potential and tunneling by approxi- 
mately I\ — 1 f=a 36%. However, the difference between the 
DFA and the true solution is not a simple rescaling of the 
axes to the nearest-neighbor interaction, as shown by the 
region bordered with dotted lines, which is the borders 
of the CP in the rescaled parameters p = p/{UI\) and 
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FIG. 2: (Color online) Phase diagram of dipolar hard-core 
lattice bosons with and without confinement-modified inter- 
action. Bordered regions represent iMPS predictions of crys- 
talline phase (CP). Power law only, free space (dashed bound- 
ary); Significant correction due to lattice, exponential plus 
power law (hatched, solid boundary); These regions use the 
parameter scalings p = p/U and t = t/U. The dotted bound- 
ary region is the confinement-modified prediction plotted in 
the rescaled parameters p = p/(UIi) and t = t/(UIi). The 
effect of the confinement-modified interaction is not a simple 
rescaling of axes. 

t = t/(UI\). Similar changes in the phase diagram will 
occur for the 2D case, with shifts of about 48% for the 
parameters of Ref. [22] . For soft-core particles which also 
possess local, isotropic interactions, the modified dipolar 
interaction will be relevant to the convexity of the inter- 
action potential and hence to the formation of supersolid 
phases [55] . 

In conclusion, we have shown that the dipole-dipole 
interaction is strongly modified by quantum fluctuations 
due to imperfect localization of particles in an unequally 
confined lattice, as found in many experiments in ultra- 
cold quantum gases. We put forward a simple character- 
ization of the modified interaction as being exponential 
at moderate separations of a few lattice sites and power- 
law for large separations. Using iMPS simulations, we 
showed that the modified interaction can significantly al- 
ter the predictions of many-body systems, including the 
determination of phase diagrams. 
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Supplemental Material for "Confinement modifies interactions of ultracold dipolar gases on optical 
lattices" 

NUMERICAL PROCEDURES TO COMPUTE THE EFFECTIVE DIPOLAR INTERACTION IN THE 

LATTICE 

In this section, we discuss two numerical techniques for computing the overlap integral 

effiyi = - 2 l d *j & fwS*) % { l*p fwV) > ( 13 ) 

f iv =wt(r) Wi/ (r), (14) 
which determines the effective interaction between particles localized in lowest band Wannier states Wj (r) . The first 



method employs the convolution theorem to express the integration over the primed coordinates in Eq. ( 13 1 as a 
function of the unprimed coordinates with high accuracy. The remaining three-dimensional integral is then performed 
with standard numerical quadrature. We will call this method the real-space method. In the second method, we 
compute the interaction matrix elements of the dipole-dipole interaction potential in the basis of lowest band Bloch 
functions. The matrix elements in the Wannier basis are then obtained by quasimomcntum Fourier transform. We 
will refer to this latter method as the Bloch expansion method. Both methods exhibit steep scaling with the linear 
domain size L which prohibits studies of very large systems. In order to ensure well-converged results, we restricted 
our analysis to a separation of at most 7 lattice sites in the main text. 



In the real-space method, we begin with Eq. ( 13 ) and apply the convolution theorem to find 

>(2), 



Gffiiil =- 2 / drf^r^iM^^Wfi^) 



where J-k[g(r)] denotes the Fourier transform of the function g(r) as a function of k and likewise for the inverse 
transform J 7 ^" "''[•]. The Fourier transform of the spatial part of the dipole-dipole potential is J r k[Cg 2 '(r)/r 3 ] = 
-47rCf ) (k)/3, @] and so 

,dd;g _ 87T 



= 11121211 



J dvf hi ,(r)F-i[CW (k)7- k [/ i2i ,(r')]]. (15) 

We now consider each Cartesian dimension to be a symmetric finite interval S = [-L/2, L/2] with periodic boundary 
conditions, and discretize each interval with n g grid points. The grid spacing in the discrete Fourier conjugate domain 
is 2ir / L and the extent of the domain in Fourier space is controlled by itn g / L, the inverse real space step size. The 
transformation from a function to its discrete Fourier conjugate is performed by the fast Fourier transform (FFT) 
algorithm in 0(n g \ogn g ) time pp. Because the Wannier functions on a finite domain are periodic and band- limited, 
their discrete and continuous Fourier transforms are related by a scaling constant provided we sample the entire 
domain at a frequency of at least twice their largest frequency component [1:. As is known for spectral methods, 



convergence of the Fourier space calculation in Eq. (15) is exponential in L provided that n g /L is large enough to 
capture the full support of the function in Fourier space. Defining g to be the support of the lowest band Wannier 
function in the discrete Fourier space, the choice n g — 2gL + 1 ensures that the function is fully captured in Fourier 
space. The support of a Wannier function in Fourier space can be determined by using Parseval's theorem on finite 
Fourier subintervals to determine that the norm is unity to a desired tolerance. For typical g ~ 5 — 7, the real space 



integration in Eq. (15) is of acceptable precision using a high-order Simpson integrator [3J. 

There are two dominant sources of error in the real space method. The first error is due to the discretization of the 
real space domain and the associated discretization error of the numerical quadrature. This error may be controlled 
by increasing g. The second source of error is spurious interactions due to periodic boundary conditions. While these 
interactions vanish as the domain becomes infinite, convergence is slow due to the power-law decay of the dipole-dipole 
interaction at long range. Hence, we have instead used a finite-size scaling analysis to extrapolate our results to the 
limit of an infinite lattice. The main limitation on the system sizes that we can reach using the real space method is 
a memory requirement which scales as O (g 3 L 3 ) due to the non-separability of the dipole-dipole potential. 

In the Bloch expansion method, we study the matrix elements of the dipole-dipole potential in the basis of lowest 
band Bloch functions ipq (r) , 

v zi = J dr J dr ' [ ^ (r) ^ (r ' )] * Mnt (r - r,) (r) ^ (r ' } ] • (i6) 
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The geometric integral Eq. ( 13 ) is then given by 

y iii 2 iiii - ip, e e e q,qi ' y > 

qiq2qiq2 

see Eq. (1) of the main text. Using the fact that the interaction potential depends only on the relative coordinate, we 
find 

yqaq 2 _ iRQyq2q 2 
q«q'i q 5 qi v ; 

where R is any Bravais lattice vector and Q = + qj, — qi — q2 modulo 2ir. This implies that there are only 

\3 • 



(L x L y L z ) independent non-vanishing matrix elements in Eq. (16) for a lattice with (L x L y L z ) unit cells, as opposed 
to the (L x L y L z ) 4 possible configurations of the four quasimomenta. As in the main text, we will study the case of 
the simple cubic lattice in which the Bloch functions separate along principal axes, 



The separability of the Bloch functions allows us to write the integral Eq. ( 16 ) in the form 



Kt= II f dr » f <[%,Aru)% 2 A<)TV int {r-r>)[% iu {r v )%^{r'J)]. (20) 

v=x,y,z~' ^° 

Changing integration variables to 2£„ = r u + r' v , 1r\ v = r„ — r' v with Jacobian 2 along each Cartesian direction and 
expanding the one-dimensional Bloch functions as 

1 1 

% (x) = lim Yl c \^ lVX . (21) 

we find 

V=X,V,Zp lu p2 U p' lu p> 3u /„={-l,l} / "''" 

Here, t is a finite Fourier cutoff used in numerics and we have defined 

i / . / . Ilv ~ 9lf + ?2i/ — fe/ , OQ s 

= Pl« - V\v + P 2 „ - V2v + Y ' ( 23 ) 

A,, = Vxu ~ Ply + P2u - P 2v H ^ ' ( ^ 

The integrals over are 

£ " / ^^e 2 ^=^, - Sin(2 7^ ) . (25) 

We will keep only the term proportional to the Kronecker delta, as this is the dominant contribution for large lattices. 
This approximation becomes exact in the limit of an infinite lattice, L — >• oo, as has been shown for the delta- function 
potential in Ref. [3]. We now also require that the interaction potential is invariant under inversion by any Cartesian 
coordinate. This is true for the q = component of the dipole-dipole interaction studied in the main text. With these 
two conditions, we find in the limit as L — > oo, 

^1 = ^3 II E * (Pl-Pi»+PL-P2» + ^ " gl »+ g '" " ^ «2<Kfont (ttA) , (26) 

v=x,y,z p lv p 2 vp' lt ,p' 2l/ 

where V (k) is the Fourier transform of the interaction potential and A = (A^, A y , A 2 ). 

The Bloch expansion method has the advantage of not introducing any discretization error, and also requires 
significantly less memory. However, the computational scaling of this method is O (i 9 ) as opposed to O (L 3 ) for 
the real space method. The Bloch expansion method suffers from spurious interactions due to periodic boundary 
conditions, just as in the real space method. We find that both methods agree when we extrapolate to the limit 
L — > oo, and allow us to put a conservative bound of 1% on our estimated error. 



DIPOLE-DIPOLE INTERACTION ENERGY OF THE ANISOTROPIC HARMONIC OSCILLATOR 

In this section we derive the dipolar interaction energy of two particles in the ground state of an anisotropic harmonic 
oscillator, Eq. (11) of the main text, and compare the results with the on-site interaction energy in the lattice. As in 
the main text, we use the quasi-2D geometry £ x = £ y = £± and choose the harmonic oscillator lengths to match the 
local curvature of a lattice site minimum via 

*» = " ( 27 ) 

Here, a is the lattice constant and V u = V v /Er with V v the lattice height as defined in Eq. (2) of the main text and 
Efj = h 2 ir 2 /2ma 2 the recoil energy. The ground state wave function may be written in cylindrical coordinates as 

1 z 2 P 2 

where p 2 — x 2 + y 2 . Using the convolution theorem, we can write the dimensionless dipolar interaction energy as 



* W ' Z) ^ eXP N^I' (28) 



rho _ » 

~2tt 2 



J dk(cos 2 0- l/3)n 2 (k) , (29) 



where 



n (k) = n (p, z) = exp -— ^ - (30) 



p 2 ^ 
4 

is the Fourier transform of the density and 9 is the angle between k and the z-axis. Performing the integration in 



Eq. ( 29 ) over z and <fi yields 



(31) 



with erf (x) the error function. Integrating over p yields 



7 ° ho -t vf m + vf kijh 2 - yl { e ± -p z f 2 cot_1 ww^m) ' (32) 

Inserting the definitions from the main text, 1= (£ 2 £ 2 L ) 1 / 3 /a, a = £ z /£j_, and (3 = a(l — a 2 ) -1 / 2 , we find 

j o ho = rril + p 2 - W - a2 y' cot_1 ^)] - ( 33 ) 

£^"k 6 



A comparison between the function Eq. ( 33 ) with the oscillator lengths chosen as in Eq. ( 27 1 and the result computed 
via Wannier functions as explained in the main text is given in Fig. [3j We find that the on-site interaction energy in 
the harmonic oscillator is always larger than the corresponding energy in the lattice, often by 10%. Additionally, the 
on-site interaction is large for small to moderate Vy , enforcing the hard-core constraint used in the many-body study 
in the main text. We stress that the on-site interaction is the only term for which the comparison with the harmonic 
oscillator is reasonable. For off-site terms, the fact that the oscillator eigenfunctions decay much more quickly than 
the Wannier functions implies significant qualitative and quantitative differences in their Hubbard parameters, as is 
already well-known for hopping [3]. 
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Planar confinement V± 

FIG. 3: On-site dipole-dipole interactions in a quasi-2D geometry as a function of the quasi-2D lattice height. The on-site 
dipolar interactions in the lattice, Iq (red solid line), and the dipolar interaction energy of a harmonic oscillator with the 
same local curvature as the lattice site, Iq° (green dashed line), have similar qualitative behavior with respect to confinement 
imbalance. In particular, both vanish as the confinement becomes isotropic, V± — V z — 45, and are significant for large 
confinement imbalance. 



